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QCD corrections to J/ijj and T production at hadron colliders 
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We calculate the cross section for hadroproduction of a pair of heavy quarks in a 3 5*1 color-singlet 
state at next-to-leading order in QCD. This corresponds to the leading contribution in the NRQCD 
expansion for J/tp and T production. The higher-order corrections have a large impact on the pr 
distributions, enhancing the production at high pr both at the Tevatron and at the LHC. The total 
decay rate of a 3 Si into hadrons at NLO is also computed, confirming for the first time the result 
obtained by Mackenzie and Lepage in 1981. 
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1. Charmonium and bottomonium states are cer- 
tainly among the most interesting systems to test 
our understanding of the strong interactions, both at 
the perturbative and non-perturbative level. More 
than thirty years after the discovery of the first 
charm- anticharm resonance, the J/ip, the study of 
their properties, including production and decay 
mechanisms, is still the subject of considerable in- 
terest 

From the theoretical point of view, a rigorous 
framework, based on the use of non-relativistic QCD 
(NRQCD) 0, has been introduced that allows con- 
sistent theoretical predictions to be made and to be 
systematically improved. However, despite theoret- 
ical developments and successes, not all the predic- 
tions of the NRQCD factorization approach have 
been firmly established. Recent measurements in 
e + e~ collisions have shown that production rates 
for single and double charmonium production are 
in general much larger than those predicted by 
leading order calculations Measurements at 

the Tevatron in proton-antiproton collisions are not 
fully compatible with those obtained at HERA in 
electron-proton collisions Q and in fixed-target ex- 
periments [Bj], suggesting the possibility that char- 
monium might be too light for the NRQCD factor- 
ization and/or scaling rules to work. In this con- 
text, a real challenge is offered by measurements of 
the J/tp polarization at the Tevatron. NRQCD pre- 
dicts a sizeable transverse polarization for J/tp 's at 
high-pT, in contrast with the latest data that now 
clearly indicate that J/ip's are not transversely po- 
larized Q. 

In view of such a puzzling scenario, it is worth to 
re-examine in detail the theoretical predictions and 
try to systematically improve on them. In this Letter 
we report on the calculation of the next-to-leading 



order correction to the hadroproduction of a pair of 
heavy quarks in a color-singlet 3 Si state. This is the 
leading term in the NRQCD expansion for a J/tp or 
a T and it is equivalent to the "color-singlet model" 
approximation. This model assumes that the non- 
perturbative dynamics leaves unchanged the quan- 
tum numbers of the perturbative quark-antiquark 
pair, which are the same as those of the physical 
bound state. It is curious to note that the analogous 
corrections to the hadronic decay width of a color- 
singlet 3 Si state, which involve exactly the same di- 
agrams, have been available for a long time [7J and 
play an important role in the extraction of as from 
T decays [8[ . Our calculation for the inclusive decay 
rate at NLO provides the first independent confir- 
mation of those results. 

2. According to the NRQCD factorization ap- 
proach Q, the inclusive cross section for direct pro- 
duction for a J/tp, T and their radial excitations, 
(hereafter denoted Q) in hadron-hadron collisions 
can be written as: 



a(pp —> Q + X) = ^ J dxidx 2 fi/pfj/ p 



(1) 
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FIG. 1: Representative Feynman diagrams for °Sp 
hadroproduction at LO (a), virtual (b), and real (c) con- 
tributions at NLO. Amplitudes with a light quark line 
(not shown) also contribute to the real corrections. 
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where p is either a proton or an antiproton. The 
short-distance coefficients a are calculable in per- 
turbative QCD and describe the production of a 
heavy quark-antiquark pair QQ state with quantum 
numbers n. Conversely, the (O|0~|O) are the non- 
perturbative matrix elements that describe the sub- 
sequent hadronization of the QQ pair into the physi- 
cal Q state. These matrix elements can be expanded 
in powers of the relative velocity of the heavy quarks 
in the bound state, v 2 . Since v 2 ~ 0.3(0.1) for char- 
monium (bottomonium) only a few terms need to 
be included in the sum over n to reach a given ac- 
curacy. Compared to the leading contribution, the 
3 S [ i ] {QQ), the color-octet terms, Osf 1 , 3 sf\ 3 Pj 8] ), 
are all suppressed by u 4 . However, the correspond- 
ing short distance coefficients might compensate this 
suppression, either by having leading order contri- 
butions starting at or s compared to the a s of the 
singlet and/or by kinematic enhancements due to a 
different scaling of the differential cross section with 
the pt- In particular it is easy to verify that the 
partonic differential cross sections at LO behave as 
1/Pt for CP-even color-octet states and \jp\ for 
the 3 S^\ compared to the 1/pf* scaling of the color- 
singlet. In the case of charmonium production at the 
Tevatron, the and l/p% behaviors are clearly 

seen in the data [9( and the introduction of the con- 
tribution from the octets allows an explanation of 
the observed pr spectrum On the other hand, 
at NLO 3 S'[ 1 ' production has both \ jp% and l/px 
components, the latter coming from the associated 
production with a heavy-quark pair, that can af- 
fect the interpretation of these behaviors in terms 
of color-octet matrix elements fllj . 

3. Analogously to orthopositronium in QED, for 
a 3 Si state the first non-zero contribution in QCD 
comes from the coupling to three gauge vectors, so 
that in hadron collisions the quarkonium is produced 
in association with a gluon, Fig. [2(a). Since we 
want to retain spin-correlations and work with he- 
licity amplitudes, we find it convenient to include 
the decay of the quarkonium state in to a pair of 
leptons, it. In so doing the projection of the quan- 
tum numbers (color and spin) of the heavy-quark 

pair to a 3 S'[ 1 ' state is automatic. The color-ordered 
helicity amplitudes for the short distance process 
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2+ 3+;5+ 6y)=0 ; 
A(l+ 2+ > 3-;5+ 67) = 
(35) 2 [12] 2 [56] 



where M = 8V2g 3 e 2 mQ. 

Next-to-leading order corrections include virtual 
corrections to the 2-^2 process, Fig. QJb), the real 
processes, such as the purely gluonic contribution 
99 — > 0,99 '1 Fig-[H C )- I n addition, channels with a 
quark-antiquark pair: gq — > Qgq and their crossings 
also contribute. 

The amplitudes gggg,ggqq — >3 S± — ► it which 
enter the real corrections, have been calculated an- 
alytically and checked numerically against those ob- 
tained with MadGraph [l2| • We note that by squar- 
ing the full amplitudes involving four gluons and 
summing over color we find a non-trivial simplifi- 
cation, with only one kinematical combination. The 
overall color coefficient is given by N C C2(F), with 
02(F) being the second Casimir of SU (3) and cor- 
responding to the color factor of the Born squared 
amplitude, i.e., C 2 (F) = [d abc / ^/W c f . 

Virtual corrections to the leading order helicity 
amplitudes have been calculated within the four di- 
mensional helicity scheme. Dimensional rcgulariza- 
tion has been adopted for isolating the ultravio- 
let, infrared and collincar singularities. The heavy- 
quark mass and wave function have been renormal- 
ized on-shell. MS renormalization has been em- 
ployed for the gluon wave function and as renor- 
malization. The color projection onto a singlet state 
leads to remarkable simplifications at NLO. We split 
diagrams containing the four-gluon vertex into a sum 
of color ordered terms, where each term is written 
as a product of color matrices in the fundamental 
representation, A a . In this way each virtual diagram 
contributes to all of the six gluon permutations each 
multiplying Tr(A a A b A c ). In fact, out of the above six 
permutations, only one independent combination is 
left which is the sum of three permutations with the 
same order of the color generators. The full am- 
plitudes are proportional to d abc / \/Nc, exactly as 
for the Born amplitude. Thus all the contributions 
from individual diagrams for each helicity configura- 
tion can be summed together, resulting in a sizeable 
simplification of analytic expressions for the final re- 
sults. We also find that the infrared divergences are 
proportional to N c , while the subleading terms in N c 
are finite, in agreement with the simplifications men- 
tioned above for the real corrections. As a further 
check of our calculation, we verified numerically the 
identity of the two color ordered subamplitudes and 
their gauge invariancc under different choices for the 
polarization vectors of the three gluons. 

To extract the singularities of the real part of the 
NLO corrections cr rcal , the dipole subtraction for- 



(Sl2 + Si 3 )(si 3 + S 23 )(si2 + S23) 



(2) 



malism 13( has been adopted, following the imple- 
mentation of MCFM [3]. The algorithm is based 



3 



LHC Vs=14 TeV 

Tevatron Vs=1.96 TeV 

p'p" - 3 s[ 1] (bb) Y(1S) + X- 

Mo-V (2m„) z + p, 2 




10000 
5000 



1000 
500 







LHC Vs=14 TeV 








Tevatron Vs=l. 


36 TeV 






p'p'- 8 sS«M 


- J/f + X — 






Ai„=V (2m,) ! + 


Pt 3 








, 1 , 


i, ,„i,„; 


1 





2.0 



FIG. 2: Variation of the cross section for direct T pro- 
duction at the Tevatron (lower curves) and the LHC (up- 
per curves) at LO (dashes) and NLO (solid). Renormal- 
ization and factorization scales are set equal, fiF = fiR, 
and no = \ / (2m b ) 2 +p?- 



FIG. 3: Variation of the cross section for direct J/tp pro- 
duction at the Tevatron (lower curves) and the LHC (up- 
per curves) at LO (dashes) and NLO (solid). Renormal- 
ization and factorization scales are set equal /j,f = Hr, 
and Ho = i/(2m c ) 2 



on mapping the singularities of the real cross sec- 
tion cr real onto an auxiliary cross section cr sub which 
can be built out of universal building blocks, the so- 
called dipoles, and it is simple enough that the singu- 
lar regions in phase space can be integrated out an- 
alytically. The remaining piece, cr real — cr sub is finite 
and can be safely integrated with standard Monte 
Carlo techniques. The divergences in the auxiliary 
cross section <7 sub are in part canceled by the the soft 
and collinear singularities of the virtual corrections 
while the rest are absorbed in the renormalization 
of the parton densities. Note that no dipoles for the 
heavy quarks are needed since there are no infrared 
divergences associated with a collinear heavy quark- 
antiquark pair in a color-singlet state. Indeed we 
find that the structure of the infrared divergences is 
the same as that of Higgs hadroproduction in asso- 
ciation with a jet [TBI, as it should be. 

4. We now turn to the presentation of the main 
features of the NLO results. In our numerical 
studies, the CTEQ6L1 (CTEQ6M) parton distribu- 
tion functions and the corresponding fitted value 
for a s (M z ) = 0.130(0.118) are used for the LO 
(NLO) predictions. The heavy-quark mass and the 
color-singlet non-perturbative matrix clement are 
(Oj^Sx)) = 9.28 GeV 3 and m b = 4.75 GcV for 
the T, and (O^ S x )) = 1.16 GeV 3 and m c = 1.5 
GeV for the J/ip. 

The inclusion of the NLO corrections leads to im- 
portant modifications of the total cross sections and 
the distributions. First we note that, even though 
the calculation of the total cross section at NLO is 
infrarcd-safc, there are some regions of the phase 
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FIG. 4: Differential cross sections for direct T pro- 

3 r ii 

duction via a S[ intermediate state, at the Tevatron 
(lower histograms) and LHC (upper histograms), at LO 
(dashes) and NLO (solid), > 3 GeV and \y r \ < 3. 
Details on the input parameters are given in the text. 



space, notably at low px and at large rapidity of 
the quarkonium state, where the perturbative ex- 
pansion is not under proper control and a fixed-order 
approximation fails to provide a reliable estimate. 
This is the same behavior found in the correspond- 
ing regions of the phase space in the calculation for 
7P -» 3 Sf ] + X at NLO A careful treatment 

of the effects of soft radiation as outlined, e.g., in 
Ref. [It} is beyond the scope of this Letter and is 
left for future work. For the sake of illustration, 
we restrict our results to the domain pj: > 3 GeV 
and < 3 where the cross section is always pos- 
itive definite. We then study how the calculation 
behaves by looking at the dependence of the cross 
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FIG. 5: Differential cross sections for direct J/ip pro- 
duction via a SJ intermediate state, at the Tevatron 
(lower histograms) and LHC (upper histograms), at LO 

(dashes) and NLO (solid). > 3 GeV and \y J ^\ < 3. 
Details on the input parameters are given in the text. 

section on variations of the renormalization and fac- 
torization scales between /j,q/2 < fiF = ^R < 2/xo, 
with /j,Q = \J (2toq) 2 + p T . The results are shown in 
Figs. Hand [3 We find that the NLO prediction for 
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r NL0 (T^Lff) =r L0 



with f3 = ±±C A - |n/ and C F = 4/3, C A = 3. The 
numerical coefficients, obtained by a Monte Carlo 
integration method, have an uncertainty better than 
one per mil, and are in good agreement with those 
of Ref. 0- Eq. provides the first independent 
check of the results of that seminal work. 

5. To summarize, we have presented the NLO cross 
section for the production of a heavy quark pair in 
a 3 S'[ 1 ' state at hadron colliders. This is the lead- 
ing contribution in the relative velocity expansion 
to the production of J/tp, T and their radial exci- 
tations. The NLO corrections are large and modify 
the expectations for the total cross sections and the 
distribution in the transverse momentum. Further 
analysis, including combination with the results of 
Ref. [Hj for pp — ► QQQ and Ref. [r|, to compare 
with the available data from the Tevatron and fixed- 
target experiments is the subject of a forthcoming 
paper. 
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T production is very well behaved both at Tevatron 
and the LHC and the scale dependence is (slightly) 
improved. For J/tp the scale dependence at NLO is 
not improved in this region of the phase space. 

fn Figs. [4] and [5] we plot the inclusive px distri- 
butions of the quarkonium states, with p^ > 3 GeV 
and |y®| < 3. The factorization and renormalization 
scales are chosen as fip = Hr = y (2t£q) 2 + pj.. 
As expected, the NLO curves have a much milder 
drop than the LO ones, due to different px scaling 
of the real contributions. This feature leads to much 
larger cross sections with respect to those predicted 
at LO, reaching one order of magnitude increase for 
p2>4(2m Q ). 

As an important check and application of our re- 
sults we also calculated the NLO corrections to the 
decay width of a 3 S'| 1 ' into hadrons. Differences with 
the production cross sections are minor and simply 
amount to the analytic continuation of the loop in- 
tegral to the physical region relevant to the decay 
process and the integration over 3-body and 4-body 
phase spaces for the virtual and real contributions, 
respectively. Our result is 
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